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INTRODUCTION 


This  report  is  a  follow-up  to  “  Numerical  Simulation  For  The  Permeation  Of  Barrier  Materials  By  Neat  Li¬ 
quid  Droplets  ”  [1],  which  presented  a  solution  to  the  governing  differential  equation  using  finite-difference  methods 
(see  also  [2]).  The  focus  there  was  upon  the  dependence  of  permeation  rate  on  barrier  thickness. 

Theoretical  results  for  the  time  dependence  of  penetrant  fluxes  were  in  qualitative  agreement  with  experi¬ 
mental  data  obtained  at  an  Army  contract  laboratory  [3].  However,  quantitative  agreement  was  not  always  satisfac¬ 
tory.  Consequently,  a  number  of  the  simplifying  assumptions  upon  which  the  model  was  based,  are  relaxed  here  in  an 
attempt  to  identify  the  sources  of  discrepancy. 

As  described  below,  there  is  significant  improvement  in  agreement  between  theory  and  experiment,  when 
account  is  t^en  of  the  finite  rate  of  penetrant  transfer  from  the  downstream  barrier  surface  to  the  bulk  sweep  gas 
stream.  This  is  primarily  due  to  the  low  vapor  pressures  of  chemical  agent  simulants,  which  minimize  the  driving 
forces  for  diffusion  through  the  gas-phase  boundary  layer  adjacent  to  the  barrier. 

Although  near-quantitative  agreement  with  experimental  data  is  achieved  in  some  cases  by  this  change  in 
boundary  condition,  there  remain  residual  discrepancies  in  the  earliest  phase  of  an  experimental  run.  Invariably,  there 
is  a  pronounced  delay  in  the  experimental  onset  of  permeation,  which  is  inconsistent  with  the  model's  predictions. 
Attempts  to  replicate  this  behavior  by  varying  the  assumed  droplet  contact  angle  and  penetrant  diffusion  coefficient, 
were  unsuccessfuL 

The  physical  basis  of  the  observed  behavior,  which  remains  unidentified,  may  also  be  the  cause  of  a  second, 
as  yet  inexplicable  observation:  the  anomalous  dependence  of  breakthrough  time  (tg)  upon  barrier  thickness  (L). 
Theory  predicts  that  tg  should  vary  as  L".  where  «  =  2.  For  some  penetrant  /  barrier  material  combinations,  n  was 
found  experimentally  to  be  as  great  as  4  or  S.  In  the  following  report,  the  attempts  to  resolve  these  issues  are 
described. 

The  model  system  considered  here  is  the  same.  At  time  t  =  0,  a  pure  droplet  is  placed  upon  an  isotropic 
membrane  in  the  form  of  a  disc  of  radius  k,  (see  figure  1).  The  contact  angle  9,  made  by  the  droplet  with  the  surface, 
is  assumed  to  remain  constant  as  sorption  proceeds  and  the  droplet  shrinks,  while  maintaining  the  shape  of  a  spherical 
section.  A  non-permeating  gas  with  zero  penetrant  concentration  sweeps  the  barrier  underneath.  The  gas  above  the 
barrier  may  also  be  flowing. 

The  dissolved  penetrant  attains  the  equilibrium  concentration,  C,  ,  at  the  droplet  base.  Whereas  the  analysis 
in  the  previous  report  was  based  on  the  assumption  of  zero  penetrant  concentration  in  the  bouom  surface,  i  =  £.,  in 
this  report  that  assumption  is  relaxed.  We  instead  examine  the  significance  of  the  finite  rate  of  mass  transfer  from 
barrier  to  sweep  gas,  and  conclude  that  it  can  indeed  be  an  important  factor  in  the  case  of  low  vapor  pressure  pene¬ 
trants. 
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Figure  1.  Schematic  diagram  of  modelled  system. 


What  follows  are  the  mathematical  fonnulation  of  the  problem,  an  examination  of  the  effects  of  finite  mass 
transfer  rates  at  the  barrier  surfaces,  concentration  dependetKe  of  the  diffusion  coefficient  in  the  membrane,  and  the 
assumed  value  of  6,  plus  comparison  with  experimental  results  attd  predictions  of  the  earlier  rruxlel.  A  preliminary 
version  of  these  results  was  presented  at  the  November  1990  CRDEC  Scientific  Conference  on  Chemical  Defense 
Research  [4]. 

Previous  work  in  this  area  also  includes  a  substantial  body  of  modelling  by  Frisch  and  coworkers  [S,  6, 7, 8, 
9],  in  which  many  of  our  observations  regarding  permeation  behavior  were  independently  made.  However,  these 
studies  did  not  focus  upon  the  issues  addressed  here. 


1.0  Mathematical  Formulation 


Following  the  earlier  model,  we  first  neglect  possible  concentration  dependence  of  the  diffusion  coefficient, 
b .  Thus,  the  concentration  of  penetrant  in  the  membrane,  d  (f,  2.  i) ,  is  governed  by  the  following  equation  in  cylin¬ 
drical  coordinates: 


dd  idC  d^C:'\ 


(1) 
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subject  to: 


C (f,  i, 0)  =  0  0<f<R,  0 <  z < £ 

(2) 

d(r,0,t)  =  t  S  0 

(3) 

=  i,,^(?.0.r)  120 

(4) 

-k,t{r,L,t)  Qir<R,  ?20 

dz 

(5) 

^  Q<i<L  r>0 

(6) 

0Szs£  ISO 

dr 

(7) 

The  caret  (^)  is  used  to  distinguish  the  above  quantities  from  the  dimensionless  ones  defined  below. 

The  km  in  equation  (4)  is  an  efTective  mass  transfer  coefficient  representing  convection  from  the  upstream 
barrio’  surface.  Conespondingly,  the  k,  in  equation  (5)  is  the  coefficient  governing  mass  transfer  from  the  down¬ 
stream  barrier  surface  to  the  sweep  gas.  The  term  A,  can  be  taken  to  denote  either  the  membrane’s  radius  or,  in  the 
case  of  a  regularly  spaced  array  of  droplets,  the  symmetry  radius  around  each  droplet  (see  figure  2). 


Figure  2.  Face  view  of  barrier  showing  a  symmetric  array  of  droplets. 

Hexagons  indicate  locus  of  symmetry  around  each  droplet 
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In  addition  to  equation  (1),  auxiliary  relationships  relate  the  time-varying  droplet  radius,  ^  (?) ,  to  droplet 
mass  losses  by  transfer  into  the  barrier  and  upstream  gas.  Thus: 


p(Vo-^0  =  + 

where: 

p  is  droplet  density  and  ('  its  volume,  which  is  related  to  A  by: 

VO)  =  ^R\})g(Q) 

and 

sin 0  (2 -I- cos 0) 

S(0)  - - - j- 

(l-(-cos0)^ 

(^0  is  the  initial  droplet  volume;  ko  is  the  initial  k  value.) 

The  cumulative  mass  flow  into  the  membrane  at  the  base  of  the  droplet,  $ ^ ,  is  given  by: 

*  o  A  f!  •  V*  .rj 

q,  =  -2Tt6jaJ/>  rz 


(8) 


(9) 


(10) 


(11) 


and  the  mass  evaporated  from  the  surface  of  the  droplet,  is  given  by: 

qE  =  k4c\\\A0)<Ci  (12) 

V 

•  P  ^ 

where  C«  =  — ,  ^  s  vapor  inessuie,  is  the  effective  mass  transfer  coefficient  for  evaporation,  and  the  exposed 

droplet  surface  area  is  given  by: 


A(r)  =  7t^^?)g(0)  (13) 

The  above  relationships  must  be  solved  simultaneously  to  determine  the  droplet  radius  and  concentrations 
inside  the  barrier  versus  time.  Once  the  concentrations  are  determined,  it  is  possible  to  calculate  the  amount  perme¬ 
ated  through  the  downstream  surface,  {qp)  and  the  amount  evaporated  from  the  unwet  portion  of  the  upstream  mem¬ 
brane  surface  (qi^)  and  the  amount  accumulated  within  the  barrier  (^^)  from: 

qp  =  2jtifc,JJ^  JJ'C  (r,  L,  I)  rdrdt  (  14) 

qM  =  {r.hrdrdi 


(15) 


The  concentration  of  penetrant  vapor  in  equilibrium  with  the  local  dissolved  concentration  in  the  top  sur¬ 
face.  is  calculated  hom: 

C  (r.t)  =  - C,  ( 1 

ti 

The  implicitly  linear  relationship  between  equilibrium  gas  and  polymer  phase  concentrations  is  an  approximation 
made  in  the  absence  of  further  data. 

In  addition,  the  amount  accumulated  within  the  barrier  is  calculated  horn: 


<lA  =  2nj^j^‘dCr,2,‘)rdrdz 


As  an  internal  check  on  the  solution,  it  must  be  true  that: 


The  solution  is  outlined  in  the  appendix.  Significantly,  the  behavior  of  the  system  is  governed  by  the  follow¬ 
ing  set  of  dimensionless  parameters: 


P 


X  =  L/Rq 


Ro 


k 

D 


kit.L 


- 

k„t.L 
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2.0  Results. 


To  validate  our  numerical  model,  we  compared  its  predictions  with  experimental  data  obtained  several  ye^ 
ago  at  Southern  Re’vaich  Institute,  in  particular  that  obtained  using  two  membrane  materials.  Neoprene  and  natural 
rubber  [3],  A  S  til  droplet  of  diisopropyl  methyl  phosphonate  (DIMP)  was  deposited  onto  10  cm^  membraites  of  var¬ 
ious  thicknesses.  In  most  experiments,  both  surfaces  of  the  membrane  were  exposed  to  an  air  flow  of  1  liter/min. 
However,  in  some  cases  there  was  no  air  flow  above  the  barrier. 

As  in  the  inevious  report  [1],  a  contact  angle  of  60P  was  assumed  in  the  absence  of  a  measured  value,  which 
leads  to  estimates  of  0.171  cm  for  Po  and  of  10.44  for  R,.  On  the  basis  of  separate  immersion  experiments  [10],  the 
Jiffusion  coefficients,  D,  of  DIMP  in  Neoprene  and  natural  rubber  were  estimated  at  7.6  xlO'*  cm^/sec  and  7.8x10'* 
cm^/sec,  respectively.  Furthermore,  based  on  the  droplet  density,  p,  of  0.98  g/cm^  and  measured  solubilities,  C,,  we 
calculated  the  partition  coefficient,  a,  (the  ratio  of  to  p)  to  equal  0.43  for  Neoprene  and  0.20  for  natural  rubber. 

The  dimensionless  mass  transfer  coefficients  at  the  bottom  surface,  k,  were  estimated  (see  appendix)  at  0.94 
and  1.83  for  Neoprene  membranes  5.6x10'^  and  1.09x10'*  cm  thick,  respectively.  The  corresponding  values  for  nat¬ 
ural  rubber  were  1.80  and  3.47  for  5.23x10'^  and  1.01x10'*  cm  thicknesses,  respectively.  For  the  experiments  per¬ 
fumed  without  air  flow  above  the  barrier  -  that  is,  with  a  sealed  upper  chamber  -  parameters  k„  and  k^  were  set  at 
zero,  since  the  amount  of  DIMP  vapor  required  to  saturate  the  chamber  represented  a  negligible  fraction  of  initial 
droplet  mass  (see  appendix).  Based  on  these  conditions,  and  the  above  (R,,  6,  a)  values,  the  theoretical  curves  in  fig¬ 
ures  3  through  6  were  generated  for  comparison  with  experimenL 


Figure  3.  DIMP  permeation  rate  vs.  time  for  a  Neoprene  membrane.  L  =  5.23  x  10'^  cm;  k^,  =  k<j  =  0 
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Figure  6.  DIMP  permeaiion  rauj  vs.  time  for  a  natural  rubber  membrane.  L  =  1.01  x  10’'  cm;  *  0 

For  both  thicknesses  and  both  materials  the  computer-  generated  curves  do  not  show  the  pronounced  delay 
observed  at  the  beginning  of  the  experimental  curves.  In  addition,  calculated  fluxes  do  not  exhibit  the  maxima  mea¬ 
sured  in  the  cases  of  thinner  barriers  (see  figure  4).  Included  in  each  figure  are  curves  generated  from  the  earlier 
model  [1],  marked  =  oo”,  which  neglected  the  downstream  mass  transfer  resistance  of  the  sweep-gas  boundary 
layer  and,  accordingly,  set  C  =  0  at  z  =  £,.  The  results  in  figures  3-6  suggest  that  fluxes  were,  in  fact,  limited  some¬ 
what  by  this  resistance. 

We  attempted  to  improve  the  fit  to  the  Neoprene  data  by  varying  O  by  a  factor  of  three  higher  and  lower  (see 
figures  7  and  8;  note:  this  caused  k,  to  vary  inversely  by  the  same  factor).  As  expected,  permeation  rates  increased  as 
b  increased  and  initial  lag  time  shortened  as  b  increased.  Significantly,  decreasing  D  improved  substantially  the  fit 
to  the  early  and  long-time  flux  data  for  both  thicknesses.  Nontheless.  manipulation  of  b  alone  is  insufficient  to  repli¬ 
cate  the  shape  of  the  flux  curve,  including  a  maximum,  for  both  cases.  Furthermore,  it  is  apparent  in  figures  4  and  6 
that  similar  adjustements  in  b  cannot  substantially  improve  the  fit  to  the  data  for  natural  rubber. 
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Figure  7.  DIMP  permeation  rate  vs.  time  for  Neoprene  membrane.  L=  5.6x10'^  cm;  kn,=kd=0;  sensitivity 
to  assumed  value  of  b. 


Figure  8.  DIMP  permeation  rate  vs.  time  for  Neoprene  membrane.  L=  1.09  xlO  '  cm;  k^=k^=0;sensiuvity 
to  assumed  value  of  b . 


Next,  the  diffusion  coefficient  was  held  constant  while  varying  the  downstream  mass  transfer  coefficient,  k,, 
by  the  same  factors  (thereby  multiplying  k,  by  the  same  factors  as  well;  see  Sgures  9  and  10).  Again  as  expected, 
increasing  the  mass  transfer  coefficient  increases  the  flux.  However,  varying  k,  from  its  estimated  value  has  no  sub¬ 
stantial  effect  on  the  permeation  time  lag.  On  the  other  hand,  decreasing  it  improves  the  fit  between  theory  and  exper¬ 
iment,  for  the  long-time  fluxes.  It  appears  that  manipulation  of  b  and/or  k,  alone  cannot  yield  quantitative 
agreement  with  the  data  for  both  thicknesses  of  Neoprene. 
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Figure  9.  DIMP  penneation  rate  vs.  time  for  Neoprene  membrane.  L-  5.6  xlO'^  cm;  k^=kd=0;sensiuvity 
to  assumed  value  of  k,. 


Figure  10.  DIMP  permeation  rate  vs.  time  for  Neoprene  membrane.  L-  1.09  x  lO'*  cm;  k„  =  =0; 

sensitivity  to  assumed  value  of  k,. 
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Finally,  we  modelled  cases  in  which  evaporation  from  the  droplet  and  upstream  barrier  surface  is  not  negli¬ 
gible,  i.e.,  where  air  flows  at  1  L/min  through  the  upper  chamber.  Values  for  the  mass  transfer  coefficient,  k„,  esti¬ 
mated  as  described  in  the  appendix,  were  2.52  for  1.02  x  lO'*  cm  thick  natural  rubber  and  1.34  for  1.09  x  10'*  cm 
thick  Neoprene.  Then,  with  the  previously  estimated  parameters  and  the  estimated  k„  and  k^  values,  the  curves 
shown  in  figures  11  and  12  were  obtained.  Comparison  of  the  experimental  curves  in  figures  6  and  12  reveals  the  pro¬ 
nounced  effect  of  upstream  evaporation  losses  in  the  case  of  natural  rubber.  (The  effect  in  the  case  of  Neoprene  is  not 
as  easily  identified  by  comparison  of  the  experimental  curves  in  figures  3  and  11,  but  is  of  similar  magnitude.)  The 
theoretically  calculated  effects  of  evaporation  from  the  droplet  and  unwet  portion  of  the  upstream  surface  are  small 
compared  to  the  effect  of  downstream  gas-phase  mass  transfer  resistance.  Much  higher  evaporation  rates  are  neces¬ 
sary  to  conform  theory  to  experiment. 


0  300  600  900  1200  1500 

TIME  (nin) 


Figure  11.  DIMP  permeation  rate  vs.  time  for  Neoprene  membrane.  L=1.09  x  lO  ’  cm.  Experiment 
with  1  L/min  air  flow  in  both  the  upper  and  lower  chambers  of  the  test  cell.  Sensitivity  to 
assumed  gas-phase  mass  transfer  coefficients. 


a  400  800  1200  1600 
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Figure  12.  DIMP  permeation  rate  vs.  time  for  Natural  rubber  membrane.  £.=1.02  x  lO  '  cm.  Experiment 

with  1  L/min  air  flow  in  both  the  upper  and  lower  chambers  of  the  test  cell.  Sensitivity  to  assumed 
gas-phase  mass  transfer  coefficients. 


3.0  Effect  of  the  Contact  Angle  (9) 

We  examined  the  sensitivity  to  contact  angle,  of  the  permeation  rate  of  DIMP  through  a  5.61  x  10'^  cm  Neo¬ 
prene  membrane.  To  do  so,  we  replaced  the  .60°  value  of  9  in  equation  (10)  with  respective  values  of  30°  and  90°. 
This  affected  not  only  the  initial  droplet  radius  (fio).  but  also  the  wetted  area  throughout  a  simulated  run.  The  result¬ 
ing  values  for  ko,  and  X,  respectively,  were  0.226  cm,  7.88  and  0.25  when  0  was  30°;  and  0.13  cm.  13.73  and  0.43 
when  0  was  90°.  Calculated  permeation  rate  curves  for  the  three  contact  angles  and  the  experimental  curve  are  shown 
in  figure  13.  Permeation  accelerates  as  0  decreases,  because  of  the  correspondingly  greater  wetted  areas.  No 
assumed  angle  yields  a  good  overall  fit.  Variation  of  0  had,  understandably,  no  effect  on  the  delay  in  the  onset  of  per¬ 
meation. 
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Figure  13.  DIMP  permeation  rate  vs.  time  for  a  Neoprene  membrane:  S.61  x  10'^  cm;  =  0; 

k,  =0.94;  a  =  0.43.  Numerical  label  denotes  value  of  8. 


4.0  Breakthrough  Time  Estimation 


In  order  to  further  investigate  the  early  time  behavior,  we  developed  an  analysis  in  which  allowance  was 
made  for  dependence  of  the  diffusion  coefficient  on  concentration  of  penetrant.  This  was  undertaken  based  on  the 
presumption  that  a  substantially  lower  diffusion  coefficient  near  z  =  £.  at  the  start  of  an  experiment,  might  explain 
the  consistently  observed  lag  in  the  onset  of  permeation. 

It  had  been  concluded  in  the  earlier  report  [  1  ],  that  theoretically  calculated  early  time  permeation  behavior  is 
frequently  indistinguishable  from  that  with  a  fully  wetted  surface.  This  observation  allowed  us  to  explore  the  ramifi¬ 
cations  of  a  variable  diffusion  coefficient  (which  complicates  the  mathematics)  in  the  context  of  a  single  spatial 
dimension  (which  requires  much  less  computer  lime  than  the  2-dimensional  model  deployed  until  now).  The  surface 
area  used  to  calculate  amount  permeated  was  that  of  the  initial  droplet/barrier  interface. 

We  again  assume  a  barrier  of  thickness  L  whose  lower  surface  is  swept  by  a  gas  with  zero  bulk  penetrant 
concentration.  However,  the  upper  surface  at  2  =  0  is  now  completely  covered  with  penetrant  and  remains  so 
throughout  an  experiment  (see  figure  14). 


13 


Downstream  Sweep  Gas 

Figure  14.  Schematic  representation  of  model  with  fully  wetted  surface. 


U 


Figure  13  depicts  typical  results  for  dimensionless  concentration  versus  dimensionless  position.  The  param¬ 
eters  correspond  to  a  fully  wetted  membrane  with  2.54  x\0'^  cm  (10  mils)  thick,  and  with  a  1  liter/min  sweep  gas 
flow  in  the  lower  surface  of  the  membiane.  The  results  from  the  two  numerical  solutions  (steady-state  and  transient) 
overlap,  as  shown  for  both  m=l  and  m=10.  Furthermore,  because  of  the  comparatively  low  estimated  value  of  k, 
(when  msl),  external  mass  transfer  (from  the  downstream  surface  to  the  bulk  sweep  gas)  is  permeation-rate-limiting, 
as  indicated  by  the  high  steady-state  dimensionless  concentration,  C,  at  z=  1  (i.e.,  most  of  the  overall  chemical  poten¬ 
tial  driving  force  is  dissipated  in  the  gas,  not  the  membrane  phase). 


Figure  15.  Steady-state  dimensionless  concentration  profiles  with  concentration  dependent  diffusion  coeffi¬ 
cient.  Diamonds  represent  steady-state  solution;  m=l.  Rectangles  represent  transient  solution  at  long 
time;  m=l.  Triangles  represent  steady-state  solution;  m=I0.  Ovals  represent  transient  solution 
at  long  time;  m=10.  L  =  2.54  x  10'^  cm.  Note  that  k,  =  0.409  when  m=l;  and  3.32  x  ICP  when 
m  =  10  (where  k,  is  defined  as  in  Eq.(23),  with  bo  replacing  b  \  the  dimensionless  analogue  of 
fcq.(29)  is  -DdCldt=k^C). 
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Having  confinned  the  validity  of  our  analysis,  we  proceeded  to  simulate  the  early  behavior  of  permeation  of 
DIMP  in  Neoprene  membranes,  as  was  first  attempted  -  and  described  earlier  in  this  report  -  using  the  2-dimensional 

model  with  constant  D .  In  particular,  we  were  seeking  an  explartation  for  experimental  results  indicating  an  anoma¬ 
lous  dependence  upon  barrier  thickness,  L ,  of  the  breakthrough  time,  /• .  defined  by  the  cumulative  permeation  of 

2 

S40  ng/cm'^.  We  now  relax  the  constant  D  assumption  by  letting  m  vary  between  zero  and  ten,  and  setting  Do  =  — 

e" 

(see  equation  29).  where  t>i  is  fixed  at  the  value  of  the  diffusion  coefficient  that  has  been  used  until  now,  which  had 
been  obtained  from  a  separate  immersion  experiment  [10].  Thus,  equation  (30)  results  in  there  being  a  lower  diffu¬ 
sion  coefficient,  at  any  concentration,  as  m  is  increased. 

We  see  (in  figure  16)  that,  when  m  =  10,  at  breakthrough  -  as  compared  to  steady-state  (figure  15)  -  the 
dimensionless  downstream  concentradon  (C  at  z=l)  is  much  lower  (0.054).  The  corresponding  steady-state  value  is 
0.644.  Thus,  the  effective  diffusion  coefficient  at  the  downstream  boundary  is  5.90  x  10'*^  cm^/sec  at  breakthrough 
and  2.16  x  10'^  cm^/sec  at  steady-state.  For  the  same  value  of  m  (10).  when  the  barrier  thickness  was  increased  from 
234  xlO'^  to  7.62  X  10*^  cm  (10  to  30  mils),  the  concentration  at  z=l  at  breakthrough,  decreased  from  0.054  to 
0.0062.  Thus,  when  L  was  tripled,  t)  (at  zsl,  at  breakthrough)  decreased  only  from  5.90  x  10'*^  cm^/sec  to  3.67  x 
10*’^  cm^/sec.  This  was  the  qualitative  effect  anticipated  when  it  was  decided  to  introduce  the  concentration  depen¬ 
dence  of  bias  L  increases,  the  effective  b  decreases,  enhaiKing  the  sensitivity  of  breakthrough  time  to  L .  However, 
because  C  (at  z^l ,  at  breakthrough)  remained  in  the  vicinity  of  zero,  the  quantitative  effect  was  marginal.  With  lower 

values  of  m,  effects  are  even  smaller.  Nonetheless,  as  described  below,  we  examined  the  theoretical  dependence  of 
on  L. 


Figure  16.  Breakthrough  time  dimensionless  concentration  profile;  m=10.  Parameters 
based  on  DIMP/Neoprene;  L  =  2.54  x  10'^  cm;  k,  =  3.32  x  10^. 


A  primary  goal  of  this  project  remains  the  prediction  of  the  relationship  between  breakthrough  time.  ,  and 

L.  Thus,  an  attempt  was  made  to  rationalize  experimental  data  which  had  previously  been  shown  [13]  to  be  express¬ 
ible  by; 

t,  =  KL  (31) 

where  k  and  n  are  constants  for  a  given  barrier^tenetrant  system.  In  one  set  of  computer  runs,  the  parameters  applied 
previously  to  simulate  DIMP  permeation  in  Neoprene  with  no  upstream  airflow  (see  section  2)  were  employed  along 

with  L  values  of  10  -  30  mils  (2.54  x  10'^  -  7.62  x  10'^  cm),  and  i,  of  either  4.52  x  10'^  or  4.52  x  10‘*  cm/sec.  The 

first  k,  value  is  an  estimated  mass  transfer  coefficient;  the  second  was  chosen  to  examine  the  effect  of  increased  mass 
transfer  resistance. 

Table  1  lists  n  values  -  obtained  from  least  squares  fits  to  the  theoretically  calculated  vs.  L  data  -  as  they 

varied  with  m,  with  k,  fixed  at  the  higher  (estimated)  value  above.  (Equation  (3 1 )  did  indeed  provide  a  good  fit.)  The 
theoretkai  results  are  in  striking  contrast  to  those  derived  from  the  experimental  data  for  DIMP  (Table  2),  which 
include  n  values  ranging  ftom  1.6  to  5.4  (1.8  to  5.4  with  no  upstream  air  flow;  1.6  to  4.2  with  an  air  flow  of  1  liter/ 

min),  for  various  barrier  materials.  Tl<e  results  obtained  using  the  k,  value  an  order  of  magnitude  lower  were  similar 

in  that  n  never  exceeded  2.  Thus,  the  concentration  dependence  of  b  cannot  explain  the  anoumalous  dependence  of 
igoaL.  Interestingly,  the  experimental  results  for  Neoprene  with  no  air  flow  above  (n=2. 14)  are  close  to  the  range  of 
theoretical  prediction.  However,  the  much  higher  n  values  for  some  of  the  remaining  materials  remain  an  enigma. 


Fable  1:  Results  of  least  squares  fit  of  n  values  (eq.  31)  to  breakthrough  times  calculated 
with  different  m  values  (eq.  30);  based  on  estimated  parameters  for  DIMP  in 
Neoprene. 

m 

n 

0 

1.66 

1 

1.74 

2 

1.79 

3 

1.85 

4 

1.90 

5 

1.95 

6 

7 

1.97 

1.985 

• 

8 

1.99 

9 

1.996 

10 

2.00 
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Table  2:  Summary  of  n  values  derived  from  least  squares  fit  of  eq.  31  [13  ]  to 
breakthrough  time  data  for  agent  simulant  DIMP  [3] 

^ith  no  upper  chamber  air  flow 

With  upper  chamber  air  flow  of  11/min 

Smithers  Rubber 

n 

n 

Butyl  ()(X)1 

No  permeation 

No  permeation 

Neoprene  (XX)5 

2.14 

1.69 

Hydrin  0(X)8 

5.40 

4.19 

SBR  0011 

2.29 

>1.53 

Natural  rubber  0010 

1.79 

1.60 

Vamac  0007 

3.03 

>1.91 

Nitrile  0004 

No  test  performed 

2.54 

Silicone  0003 

No  test  performed 

2.31 

5.0  Conclusions. 


Improvement  in  the  fit  of  modelling  results  to  experimental  data,  has  been  achieved  by  inclusion  of  down¬ 
stream  gas-phase  mass  transfer  effects.  However,  there  remain  marked  discrepancies  between  theory  and  observation 
in  the  case  of  early-time  permeation  behavior,  leading  up  to  the  “breakthrough  ume”: 

i)  The  experimentally  observed,  pronounced  delay  in  the  onset  of  permeation  remains  irreconcilable  with  the 
model,  even  after  including  the  downstream  gas-phase  mass  transfer  resistance,  as  well  as  a  concentration-dependent 
penetrant  diffusion  coefficient  in  the  barrier,  and  adjustment  of  droplet  contact  angle. 

ii)  Similarly,  variation  of  model  parameters  -  in  particular,  those  governing  concentration  dependence  of  the 
diffusion  coefficient  -  proved  unsuccessful  in  replicating  the  experimentally  observed  variety  of  dependences  of 
breakthrough  time  upon  barrier  thickness. 

This  leads  us  to  conclude  that  either  the  experimental  data  -  at  least  at  early  times  -  were  not  accurately  mea¬ 
sured;  or  physical  phenomena  -  e.g.,  complications  arising  from  the  presence  of  chemical  additives  in  as-received 
rubber  samples,  or  a  non-equilibrium  time -dependent  droplet  contact  angle  (droplet  spread)  -  are  responsible  for  the 
observed  dynamics  of  permeation.  The  immediate  plan  is  to  attempt  to  reproduce  the  anomalous  experimental  results 
for  at  least  one  penetrani/barrier  material  combination. 
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I  List  of  Notations _ 

A  Droplet  exposed  surface  area. 

to  Equilibrium  concentration  at  the  droplet  base. 

C  Dimensionless  concentration  as  defined  in  the  appendix  by  equation  (34) 

.V 

C  V)qx)r  concentration. 

.V 

Ct  Equilibrium  vapor  concentration. 

(f  Finite-difference  approximation  for  the  intermediate  value  which  arises  from  the  implicit  computation  of  C. 
b  Solvent  diffusion  coefficient 

Gas  phase  diffusion  coefficient. 
k  Proportionality  constant. 

kd  Mass  transfer  coefficient  representing  evaporation  from  the  surface  of  the  droplet. 
k^  Dimensionless  mass  transfer  coefficient  as  defined  by  equation  (39). 
k„  Mass  transfer  coefficient  representing  convection  from  upstream  barrier  surface. 
k^  Dimensionless  mass  transfer  coefficient  as  defined  by  equation  (40). 

k.  Mass  transfercoefficient  governing  transfer  of  mass  from  downstream  barrier  surface  to  the  sweep  gas. 
k.  Dimensionless  mass  transfer  coefficient  as  defined  by  equation  (38). 

L  Barrier  thickness. 

nid  Initial  droplet  mass. 

P  Partial  pressure. 

V 

P  vapor  pressure  of  the  solvent  penetrant  which  was  diisopropyl  methyl  phosphonate. 

P,  Total  pressure. 

9a  Amount  accumulated  within  the  barrier. 

Dimensionless  amount  accumulated  within  the  barrier  defined,  in  equation  (44). 
qa  Cumulative  mass  flow  at  the  base  of  the  dropleu 

qg  Dimensionless  mass  flow  at  the  base  of  the  droplet  as  defined  by  equation  (42). 

9£  Mass  evaporated  fiom  surface  of  the  droplet 

qg  Dimensionless  mass  evaporated  from  surface  of  the  droplet  as  defined  by  equation  (43). 

9j*f  Amount  lost  from  unwet  portion  of  the  upstream  membrane  surface. 
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q£  Dimensionless  mass  evaporated  from  surface  of  the  droplet  as  defined  by  equation  (43). 
qi^  Amount  lost  from  unwet  portion  of  the  upstream  membrane  surface. 

qj^  Dimensionless  amount  lost  from  the  unwet  portion  of  the  upstream  membrane  as  defined  by  eq  (45). 

Amount  penneated  through  the  downstream  surface, 
qp  Dimensionless  amount  permeated  through  the  downstream  surface,  as  defined  by  eq  (4 1 ). 
r  Radial  coordinate, 

r  Defined  as  r/Ro 

R  ( t)  Time- varying  droplet  radius. 

Re  Reynold’s  number. 

R,  Membrane  radius. 

Rj  Dimensionless  membrane  radius  as  defined  by  eq  (32). 

Sc  Schmidt  number 

Sh  Sherwood  number. 

?  Time. 

t  Dimensionless  time, 

tp  Breakthrough  time. 

V'  Droplet  volume. 

^0  Initial  droplet  volume, 

z  Distance  from  the  upstream  barrier  surface, 

z  Define  by  equation  (35). 

a  Defined  after  equation  (76). 

^  Defined  after  equation  (76). 

e  Contact  angle. 

X  Defined  following  equation  (46). 

p  Droplet  density, 

o  Defined  by  equation  (37). 
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Appendix 


8.0  Conversion  to  Dimensionless  Variables 

To  identify  key  parameters  and  generalize  the  results,  equations  (1)  -  (18)  were  rewritten  in  terms  of  the  fol¬ 
lowing  dimensionless  variables: 


kjC  X 
DC. 


t>t. 


(32) 

(33) 

(34) 

(35) 

(36) 

(37) 

(38) 

(39) 

(40) 

(41) 
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<tB  =  — 
ntj 

(42) 

<?£ - T" 

mj 

(43) 

>> 

II 

(44) 

K 

II 

(45) 

ntj  is  the  initial  droplet  mass. 

Equation  (1)  becomes: 

1  3C  _  3^C  1 3C  1  3^C 

1}^  ~  3?  ^  ^  3z^ 

(46) 

where  X  =  L/Hq. 

The  initial  and  boundary  conditions  become; 

C(r,z.O)  =  0  0<rS^,  OSzS  1 

(47) 

C(r.0.0  =  1  OiriRU)  t>0 

(48) 

dC(r,0,i) 

=t^C(r.O,0  R{t)<r<R^  t>0 

(49) 

dC  (r,  1,  t) 

—  =  -t,C(r,  1.0  0<r<R^  ;20 

(50) 

dC(/f,.z.  0 

- if -  =  0  0<z<l  f^O 

dr 

(51) 

dC  (0,  z,  1) 

\  -  0  0<z< 1  f>0 

dr 

(52) 

Furthermore,  equations  11,12,14,15  andl7  become,  respectively: 

6oX  ri  r*(i)  3C(r,  0,  f) 

"'’-scejloJ.  a, 

(53) 

0 
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3oX)S:^ 


6<jXJk  fg 


The  integrals  were  evaluated  by  applying  Simpson’s  rule  to  the  numerical  values  of  C(r,z,t)  which  had  been  deter¬ 
mined  as  described  below. 

9.0  Calculation  of  the  Mass  Transfer  Coefficients 

In  order  to  obtain  a  representative  value  for  k,,  we  referred  to  the  configuration  of  the  liquid-droplet  chal¬ 
lenge  permeation  tests  conducted  at  Southern  Research  Institute  [3].  In  these  experiments,  a  cylindrical  test  cell  was 
divided  into  upper  and  lower  chambers  by  the  permeation  barrier.  In  the  lower  chamber,  a  Teflon  insert  was  used  to 
accelerate  the  air  stream,  thereby  promoting  gas-phase  mass  transfer  (see  figure  18).  The  volume  of  the  lower  cham¬ 
ber  was  6  cm^  without  the  Teflon,  and  3.5  cm^  with  the  Teflon  inserted.  The  volume  of  the  upper  chamber  was  16 


AGENT  OR  SIMULANT  DROP 


TEST  SAMPLE 


TO  CHARCOAL 
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SYSTEM 


TEFLON  INSERT 


FIGURE  18.  Experimental  configuration  of  multichamber  test  cell  [3] 


The  mass  transfer  coefficient  at  the  bottom  surface  was  obtained  from  the  following  correlation  (14): 
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In  addition,  Q=air  flow  rate,  A=  channel  cross-sectional  area  of  test  cell,  p=density  of  air,  d=  channel  depth 
below  sample,  |i= viscosity  of  air,  T=absolute  temperature  in  Kelvin,  R=  ideal  gas  constant,  Pi=absolute  pressure  in 
atmosphere,  P=  partial  pressure,  1cg=  gas  phase  mass  transfer  coefiBcient  and  for  simplicity,  carets  (^)  have  been 
omitted  above  the  symbols  for  dimensioned  parameters.  To  estimate  the  Reynolds  and  the  Schmidt  numbers,  we 
refer  to  the  experimental  conditions.  Air  at  25*^  C  and  1  atmosphere  was  fed  to  the  downstream  chamber  at  1000  ml/ 
min.  The  cylindrical  cell  which  contained  the  sample  had  a  diameter  of  3.57  cm.,  a  height  above  each  sample  of  1.59 
cm.,  and  a  depth  below  each  sample  of  1.27  cm.  (see  figure  19). 


FIGURE  19.  Muliichambcr  test  cell 


The  gas-phase  diffusion  coefficient  (Dy^g)  was  esfimated  using  [IS] 


P[av)Y^+  (Zv)i^’]^ 


where  A  and  B  were  air  and  diisopropyl  methyl  phosphonate  (DIMP).  respectively,  and: 
M  s  molecular  weight 


V  =  atomic  diffusion  volume 


was  thereby  estimated  to  be  6.34  x  10'^  cm7sec. 


(59)  to: 


Once  the  Sherwood  number  was  obtained,  we  derived  the  mass  transfer  coefficient  io  by  rearranging  eq 


f  _ 

RTDP 


The  dimensionless  expression  for  k,  is  given  by: 


where; 

kef’'' 

k,:  the  mass  transfer  coefficient  defined  by  eq.  (5),  is  equal  to - 

L:  thickness  of  the  membrane. 

V 

I*  :  vapor  pressure  of  the  solvent  penetrant,  diisopropyl  methyl  phosphonate  (DIMP),  0.27  mmHg. 

C,  :  equilibrium  dissolved  concentration  at  the  droplet  base. 

t):  diffusion  coefficient  of  the  solvent  in  the  barrier  material.  Its  value  was  obtained  from  separate  immersion  exper¬ 
iments  [10]. 


The  value  for  the  gas-phase  mass  transfer  coefficient  above  the  unwet  surface  of  the  membrane,  was  calcu¬ 
lated  using; 


where  the  Sherwood  number  (Sh)  and  diffusion  coefficient  (D^g)  ^  estimated  as  before. 


The  mass  oansfer  coefficient  was  made  dimensionless  as  follows: 


where: 

.V  V 

C,  is  vapor  concentration  in  equilibrium  with  the  droplet,  P  /RT. 


(66) 


In  experiments  with  no  air  flow  in  the  upper  chamber  of  the  test  cell  (see  figure  18).  the  Reynolds  number 
(Re)  was  set  to  zero  (V=0  in  equation  60),  which  gives  a  limiting  value  of  0.43  for  the  Sherwood  number  (Sh).  Sub¬ 
stitution  of  this  value  into  equation  65,  results  in  a  mass  transfer  coefficient  (k^)  of  7.64  x  10'^  cm/sec.  The  corre¬ 
sponding  dimensionless  !(,„  (equation  66)  was  sufficiently  close  to  zero  (0.0170)  to  justify  neglect  of  evaporation 
fiom  the  upstream  surface.  (Note,  in  addition,  that  P  is  so  low  as  to  justify  neglect  of  the  amount  of  mass  lost  from 
the  droplet  in  saturating  a  closed  upper  chamber. )The  same  value  (0.0170)  was  obtained  for  Ic^. 


10.0  Finite-Difference  Approximation. 


Figure  20.  Arrangement  of  grid  for  finite  difference  analysis  (note:  grid  is  finer  than  shown;  Iq  is  not  the 
grid  point  after  0) 


2.8 


The  finite-difTerence  method  [16. 17, 18, 19]  was  used  to  obtain  the  concentration  profile  C(r,z,t)  in  the 
membrane.  To  solve  the  pardai  differential  equation  (46)  governing  concentration  using  the  finite-difference  tech¬ 
nique,  the  membrane  was  first  divided  into  grid-points  (i  J).  denoting  space  points  having  coordinates  iA/- ,  yAz .  To 
minimize  computation  time,  the  implicit  alternating-direction  method  developed  by  Peaceman  and  Rachford  [  14]  was 
again  used  to  obtain  the  numerical  approximation  for  the  concentration  profile  C(r,z,t). 

The  method  consists  of  alternatively  treating  the  respective  spatial  derivatives  in  the  r  and  the  z  directions  as 
unknowns  in  successive  dimensionless  time-steps  A//2.  The  first  half  time-step  is  implicit  in  the  r-direction.  and  the 
second  half  time-step  is  implicit  in  the  z-direction  [17].  The  net  result  is  the  concenuation  C(r,z.t)  at  the  end  of  inter¬ 
val  A/.  If  we  denote  the  set  of  dimensionless  concentrations  at  “old”  time  as  Qj^,  those  at  the  end  of  the  first  half 

time-step  as  C^  -,  and  the  “new"  values  (at  the  end  of  interval  A/)  as  ,  the  finite-difference  analogs  to  equation 

(46)  become: 


~  n  _  1 2  r  Q-l./  1./  1./ ~  ^i- 1./ 1,  « 

(Ar)/2  1,  2i(Ar)^  J  (Az)^ 

^i.y.ii+t  ~  ^i,;  _  ^2  l.y  ^  ~  ^  -  --I,**  1  ~  ^i.j*  l,M*l 

(^0/2  V  (Ar)^  2i(Ar)^  J  (Az)^ 


(67) 

(68) 


Equations  (67)  and  (68),  plus  the  boundary  conditions  outlined  below,  each  fonn  a  uidiagonal  matnx  of 
equations  in  tenns  of  unknown  concentrations.  In  order  to  obtain  Table  3,  equations  67  and  68  were  used  with 
1  s  i  s  /j  -  1  and  I S;  +  1 .  Due  to  the  specisd  boundary  conditions  (sec  next  section)  that  exist  at  j=0  (z=0)  and 
j=JL  (z=l).  and  i  =  0  (nsO)  and  i  =  Ij  (rsR,)  the  tridiagonal  mauix  equations  obtain  by  using  equations  (67)  and  (68) 
were  not  directly  applicable  at  these  positions.  Separate  matrices  of  equations  were  necessary  in  order  to  apply  these 
boundary  conditions.  To  solve  the  tridiagonal  matrices  in  Tables  3  through  8,  we  used  the  Thomas  algorithm  [18]. 


10.1  Boundary  Conditions. 

At  the  unwet  portion  of  the  upper  surface,  i.e.  j  =  0,  i>Io  (see  figure  20),  boundary  condition  (49)  applies.  In 
the  preceding  report  [1],  a  “symmetry”  analog  was  used  to  eliminate  the  virtual  concentration  C_j  ^  in  the  finite  dif¬ 
ference  approximation  to  the  derivative  in  equation  (49).  However,  that  is  justifiable  only  when  (dC/dz)  is  zero.  For¬ 
tunately,  that  was  generally  true  of  the  cases  examined  in  that  report.  For  greater  generality,  we  employ  the  “quarter 
point”  approach  i.e.,  we  write  the  partial  differential  equation  at  J=l/4,  making  the  following  approximation: 


3C 

^  3  8C 

18C 

di 

i,j  »  1  /4 

0 

1./  s  0  1 

d^C 

^  3  8V 

^  1  8V 

dr^ 

1/4  ^ 

1,/  »  1/4 

4TT 

../-o 

dC 

,  38C 

^  18C 

dr 

1/4  ^ 

1,  /  =  1  /4 

*  0  1 

(69) 

(70) 

(71) 
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Az 


I  1,/  3  1/4 


where: 


^*.1  ^*.0 


I  *.7  «  1/2 


!£l  .  *.c 


(according  to  the  boundary  condition). 


We  substituted  the  above  approximations  in  equation  (46).  What  follow  are  the  general  equations  for  j=0. 


when  /q  <  i  s  /j : 


3(^-l)C;.,.o  +  6(l+p)C*o-3(l  +  i)C;,,.o 


(i-i)C,.,-2(i+P)c’,+(i  +  i)c;,,,, 


+  (8a+2p)C,,.,+  l6p-8o(l  +  i,Az)]C.o., 


[6P  +  8a  ( 1  +  k^^z)  ]  C,. 0. -  (8a  -  2p)  q 


3(1 -^)c;.,.o+6(p- DC"  0+3(1 +  i)c;,,_o 


+  (i-^.)c;.,,,+2(p-i)c’,  +  (i  +  i)c;,,  J 


1  .4  1 

where  P  =  — — —  and  a  =  ^  (— ) 
At  X^'Az' 
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A  funher  complication  arises  from  the  boundary  condition  at  r=0  (i=0).  Here  th^  indeterminate  form  0/0 

1  do  1  do  d^O 

results  for  -  -5- .  However,  by  applying  I’Hdpital's  rule,  -  3-  becomes  — - .  Therefore,  at  i=0  the  partial  difTerendal 

f  dr  r  of  0^2 

equation  becomes: 


ac  a^c 

ar^  ar" 


(77) 


and 


The  finite-difTerence  analogues  to  (77)  at  z=0  (j=0)  (see  figure  20)  are  dien: 

(6p+12)Co,o-12Ci.o  =  -2(p  +  2)Co.i+4cI_,+  (6^  -  8a(l  +  ifc^Az) )  Qo.,  +  2  (4a+ P)  Q.  1,,  (78) 

(2P-8o)Co.,..^,+  [6p  +  8a(l+i(:,Az)lCo.o..*,  =  12C;.o+  (6p- 12)C;,,  +  (2p-4)C;_,  (79) 

ac 

At  isli  (r=R^,  3-  is  again  =0.  As  result,  here  the  finite-difference  analogs  to  (46)  are; 

of 

-6C;^-i.o  +  6(l-t-p)C)^.o  =  2C;^-i.i-2(l  +  P)C;,,i+(2P-H8a)C,^_,  ,+  (6p- 8a(  1  +  ( »<>) 


and 


[8a(l  +  **Az)  +6p]C,^_o..+  i+  (2p-8a)C,^_,  = 


- 1.0  ■•■  ^  (^  ^/„o  .  I  , 


+  2(p-  DC,  , 


(81) 


As  emphasized  earlier,  a  difTerent  boundary  condition  is  applied  at  j=jL  (z=I)  than  that  used  in  the  previous 
report  (see  equation  SO).  In  order  to  incorporate  it  into  the  finite-difference  analysis,  we  again  used  the  quaner  point 
approach.  Thus,  the  second  partial  derivative  with  respect  to  z  is  approximated  as: 


where: 


ac 


ac 


(■3'>  -(3-> 

t,,L  t./L-\/2 

Az 

T 


ac 


=  -k.C. 


'./t- 


‘./t- 


(82) 


(83) 
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(according  to  boundary  condition  (SO)),  and: 


<3“)  “* - JZ - 

02  ;/  _  1  yo  ^2 


i.jL-\n 


The  finite-difference  analogues  at  j=jL  (z=l)  then  become: 


(*  “ 57) 2(l+P)C;y£_| 


3(^  -1) c\- X.JL  +  6 ( 1  +  P) - 3 ( ^.  +  1) c., -  +  ( ^.  +  1) c*, ,  +  (8a  +  2p) , 


for  the  first  half  time-stq)  and: 


[6P-8a(l+A:,A2)]C,,^i,, 


(2p-8a)q^^.,,,^,+  [6p  +  8a(l-t-*,A2)lC..^t,,^,  = 


+  3(i  +  l)C;„.^^+(l-i)C;.,.,^.,  (86) 


-K2(p-l)q^i.,+  (^,  +  l)C.,,.,^., 


for  the  second  half  time-step. 

The  convection  boundary  condition  (equation  50)  applies  over  the  entire  downstream  surface.  Accordingly, 
equations  (85)  and  (86)  were  applied  at  all  values  of  i  except  at  the  end  points  (i=0  and  i=Ii),  where  the  appropriate 
boundary  conditions  (equation  51  and  52)  were  again  applied.  The  result  is: 


(6p-H2)Co.^^-12Co.,i  = 


-2(P  +  2)Co^^_,  +4C,  -K2(4a+  PlCo^^t-i,, 

+  [6p-8a(l+l:,Az)]Co.,i,, 


"  -  l./t  ^  (  P  ^  )  ^l,.jL  ~ 


2  (P  2P)  1,, 

+  [6p-8a(l 
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for  the  first  half  time-step,  and: 


^O.yL  ■*■  ^^l.yL- 1’ 

■*■  (2P  ~  4)  Cq  j[^_  I 


(89) 


(2p-8a)Co,;i_,.,^,-K  [6P-h8a(l-t-l:,Az)]C,_.^y.  = 


r^Q,  -  l.yi,  ^  (P  '  )  ^/,,yt  -  l.jL  -  I 


+  2(P-l)Cy,_,^.. 


(90) 


for  the  second  half  time-step 


11.0  One  Dimensional  Transport  with  Concentration  dependent  Diffusion 
Coefficient. 


For  the  fully-wetted  barrier  with  concentration-dependent  D,  the  dimensionless  partial  differential  equation 
becomes: 


-  n^£  — 
dt  32^  ^  dz 


where 


Since  eq  (91)  is  non-linear,  we  linearize  it  to  become: 


dC 

It 


dD  '* 


dC 


dz 


(91) 


(92) 


where  y  denotes  conditions  at  “old  time”  -  i.e.,  at  the  start  of  a  finite-difference  time  step.  To  arrive  at  the  solution 
presented  in  Table  9,  the  finite-difference  approach  as  described  in  the  previous  section  was  used.  The  principal 
change  (other  than  in  dimension)  was  that  the  boundary  condition  at  J=jL  was  applied  using  the  following  variation  of 
the  1/4  point  approach: 


dC  dC  1  fd^C) 

(^)  -(^)  -jAz  ^ 

yt.-I/4,«  +  l  OZ  4  V3z 


(93) 
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where: 


dC  K  ^ 


m. 


ac  ac^ 

//,*+!  ;£-l/2.»-»l 


yL-l/4,«+l 


/L-l/2.<i+l 


Az 


It  follows  that: 


(a. 


2I>'' - 1 . . .  I  -  2  a, Az  +  D'')  „ 


jL~l/4,i>*l 


Table  3: 2-Dimensional  Analysis:  TruUagonal  Matrix  Equations  for  the  First  Time-Step  (implicit  in  r  only) 
with  j*Q;  j*jl  and  . 


(2P  +  4)C  0./-4C  1,; 


i  =  1  (^.-l)C*oj  +  2(l  +  P)C’i.,-  (l  +  i)C\, 


i  =  ...  (l-l)CVl.;  +  2(l  +  P)C,.,-  (1  +  ;^)C.>,., 


»  =  ^-1  (;J;-1)C’,,-2.>  +  2(1  +  P)C’,,-u-  (l  +  i)C‘/,., 


«■  =  /,  -2C/,-,.;  +  2(l  +  P)C/,.; 
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di  =  aq^_,  ,  +  2(P-a)q^  ,  +  aC,^^,,, 


ify^O 


di  =  (l-^)Ci.i.i-2(l  +  p)C\i+  (l  +  ^)Ci^i.i+  (8a  +  2p)C,.,  ,+  [6p- 8a(  1  +  i^Az)  ]  qo., 

if  i  =  0;  j>/o 

Table  4;  2-Dimtnsional  Analysis;  Tridigonal  Matrix  Equations  for  the  First  Half  Time-Step  (implicit  in  r  only) 
with  j^O;  lQ<i^  1 1 


i  =  /o+l 

6(1  +  P)C/,+i.o-3(l  +  i)C/,+2,o 

i  =  /o  +  2 

3(^-l)C /,♦  1. 0  +  6(1  + P)C/,  +  2,o-3(1  +  ^)C/,  + 3,0 

=  di^+2 

i  =  ... 

3(^.-l)C/_,.o  +  6(l  +  p)C\o-3(l  +  i)C.*,.o 

=  d, 

i  =  /,-l 

3(i.-l)C,,-2.o  +  6(l  +  p)C,.-i.o-3(l  +  ^)C/,.o 

=  d,., 

i  =  /, 

-  6C  /,  - 1. 0  +  6  ( 1  +  P)  C  0 

di  =  (l-^)C’i_i.i-2(l  +  P)C,,+  (1  +  1)C.>,.,+  (8a  +  2P)q,  ,+  [6p  -  8a  ( 1  +  ^Az)  ]  C,  „  , 


for  j  >  /q  +  1 

1)  forj  =  /o+l 

d,^  =  2C’/,  - 1. 1  -  2  ( 1  +  P)  C /,.  I  +  (2P  +  8a)  ,  +  [6P  -  8a  ( 1  +  k^Az)  ]  q_o  for  i  =  /, 

do  =  -2(P  +  2)C,,i  +  4C.  +  ,.,  +  [6P-8a(l  +k„Az)]C,o+  (8a  +  2P)q,  fori=0 


J5 


Table  5:  2-Dimensional  Analysis:  Tridiagonal  Matrix  Equations  for  the  First  Half  Time-Step  (implicit  in  r  only) 


with  jsUL:  0  5 » s  . 


i  =  0 


(6P+12)Co.yi,-12C,.ji. 


= 


*=1  +  6  ( 1  +  3)  C  \  jL  -  3  (^.  +  1)  C  2,yL  = 

i  =  ■■■  3  (^  -  1)  C  ,_i,yi+ 6(1  +  P)C  3  (^  +  1)  C  i+i.y/,  =  d^ 

i  = /[  -6C/,_i,7/,  +  6(1 +  p)C/,.yf,  ~ 


do  =  -  2  (P  +  2)  C” o.yt  - 1  +  4C  i.yz.. ,  +  2  (4a  +  P)  Co^y^.. ,  +  (6p  -  8a  ( 1  +  k,Sz)  ] 

dj  =  ~  Yi^  ^  i~2(l+P)C,.yi,-i  +  ^Yi*  +  (8a  +  2P)  C,  y^_  i  +  [6P  -  8a  ( 1  +  k^Az)  ]  C, 

d,^  =  2C  /,-i.yL-i  -  2(P  +  1)  C  /,.y/,-i  +  (8a  +  2P)  +  [6p  -  8a(  1  +  k,Az)  ]  C,  ji^ 


Table  6:  2 -Dimensional  Analysis;  Tridiagonal  Matrix  of  Equations  for  the  Second  Time-Step  (implicit  in  z) 
with  /q  +  1  5  i S  /,  if  Iq*0 


e 

J  =  0 

[6p  +  8a(l+it^Ar)}qo 

-(8o-2P)q,  =do 

j  =  1 

0,  *«+  1 

+  2(P  +  a)q,,^,  -aC„2_,^,  =di 

j  =  ••• 

l.*  +  1 

+  2  (P  +  Ct)  C,j_  J  ~®^i,>+l,/l+l  ~ 

1 

II 

+  2(P  + a)C, ~  dji^.i 

dj  =  (1  -  ^)C  i_i,y  +  2(p  -  1)C  +  (1  +  +  for  i ?t/,.  j  ^0  & 

d’j  =  2c'/,-i,;  +  2(P-  l)C’/,.>  fori  =  /,;;*0 
d f  =  (2P - 4)  C  0,/  +  4C  ij  for  i  =  0\  j*0 

d  j  =  6C  i-i,o  +  6(P  —  1)C  j,  0  +  2C  i-1,1  +  2(P  —  1)C i  for  i  =  l\  \  j  =  0 

d  j  =  12C  1,0  +  (6P  —  12)  C  0,0  +  4C  i,  i  +  (2P  —  4)  C  o,  i  for  i  =  O',  j  =  0 

dj  =  3(l-i)C’i-i,o  +  6(p-l)Ci.o  +  3(l  +  i)Ci>,.o+  d  -  +  2(P- 1)C  ,,  +  (1  +  1)C 


for  i*l^,  i*0  &  j=0 
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Table  7;  2-Dimensional  Analysis;  Tridiagonal  Matrix  Equations  for  the  Second  Time-Step  (implicit  in  z) 


with  0  S  i  S  /,  and  j  =  JL 

=  JL  (2P  -  8a)  Ci  jL +  [63  +  8o  ( 1  +  =  dji 

d  JL  ~  3(1-  yP  ^  i-\,JL+  6(3  ~  1)  c  ijL  ~  Yi^^  1  ■*■  2  (3  ~  1)  c  ,,JL-l 

+  (^  + for  iTtO  &  i5t/, 
d  JL  —  12C  1,7/:.  +  (63  —  12)  C  o.7i  +  1  +  (23  —  4)  C  o,yi- 1  fori=0 

d  JL  =  6C  /, -1.71  +  6(3-  1)C  i^,jl  +  2C  i^-i,jL-i  +  2(3-  1)C  /„7i.-i  fori=li 

Table  8:  2-Ditnensional  Analysis;  Tridiagonai  Matrix  of  Equations  for  the  Second  Time-Step  (implicit  in  z) 
for  Omi^.j^Q 

j-^  ~  ®^i.0.«+l  ^  l.n+I  “  ®^i.2,»+ I  -d] 

J  —  •••  “  *^^i,7- 1,*+ 1  ^  (P  ^i.7.»  +  l  ~  ®^i.7  + 1.(1+ 1  ~ 

~  ®Ci.7-3,  *♦  1  2  (P  +  Ol)  C,j_2,,+  i  —  OCj_y_  J  I  =  d  j-2 

j  =  JL  (23  -  8a)  c,  7t_,  +  [63  +  8a(l  +  k^Az)  ]  C,  ji^  =  d ,, 
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for  2  <  y  <  y  -  1 


d,  =  (l-^)C‘i-i.>  +  2(3-l)C’,y+  (l  +  ^)C'.*,j 
dj  =  (2P-4)C  o,y  +  4C  i.y 
dj=  (l-l)C*i-,.y  +  2(P-l)C\,+  (l  +  i)C..u+a  forj=l 

d  JL  =  3(l-^)C,_i,7i.+  6(P-l)Ci_yt  +  3(^.  +  l)C/+i_yt+  (l“^)Ci-i,yi,  +  2(P-l)C,,yt-l 


+  (1  +  1)  Ci^yjL-X  =  yi  &  i  ^  0 

djL  =  12C*,.yt+  (6p-12)C*o.yt  +  4C‘i.//.-i+  forj=JL&i=0 

Table  9:  Tridiagonal  Matrix  Equations  For  The  One-dimensional  Fully-Wetted  Surface  Model. 


j=l: 


(2X1)''+ _^(C2.,-l)  +X2?’' 


'2,11+1 


= 


j=2: 


\dD^ 


(C3.,-C,.,)-XD’' 


C,,^,  +  (1+2X£»'")C2.,,,- 


XdD^ 
L  4 


^3,«+l  •"  ^2 


! 


^/  +  1,  «  +  I 
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j=jL-l: 


XdD* 


(C>,  ,-C,, 


'-JL-2,m. 


rxdD^  ^  ^  j 


j=jL; 


1  -  8XD’'  -  2A/<iD’' 


C,/.- 1.-.  1  +  [3  +  8X(Azi,  +  D’')  ]  ,  + 


2A/dD 


,»  «)  -  2Azk^) 


(D'^^2) 


^jL,m+l  -  ^jL 


kD 


=  C,..--^(C2,.-1)+XD' 


d  =  C 
1 


^JL  -  ^;Z,-l.«+ 3Cy£^, 


<iD'''  is  the  derivative  of  o'*  with  respect  to  C. 
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12.0  Programs  Listing. 


12.1  2*  Dimensional  model 
PRCXSRAM  DRPEV 
C 

C  DROPLET  DIFFUSION  THROUGH  A  MEMBRANE  (INCLUDING  EVAPORATION) 

C  FINITE-DIFFERENCE  SOLUTION  OMPLICIT  ALTERNATING  DIRECTION  METHOD) 

C  SIMPSON’S  RULE  IS  USED  TO  INTEGRATE  OVER  AREAS 

C  TRIDAG:  SUBROUTINE  FOR  SOLVING  A  TRIDIAGONAL  SYSTEM  OF  SIMULTANEOUS 
C  EQUATIONS. 

C  CPRIME;  VECTOR  FOR  TEMPORARY  STORAGE  OF  CONCENTRATION  COMPUTED  BY  TRIDAG 

C  CSTAR:  MATRIX  OF  CONCENTRATION  C*  AT  THE  END  OF  THE  FIRST  HALF  TIME-STEP 
C  IFREQ:  NUMBER  OF  TIME-STEPS  BETWEEN  SUCCESS  IVE  PRINTING  OF  CONCENTRA- 
C  TIONS. 

C  AMS3MS,CMSJEMS  JTV1S,GMS.AM3M.CMJ>M  ARE  THE  COEFFICIENT  VECTORS. 

C 

C  IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

CHARACTER*  15  FNAME 

REAL  LAMBDA.KMJCD4CS.Q0.QACC,QT,QED 

DIMENSION  AM(0:350)3M(0:350),CM(0:350)J)M(0:350),C(0:350.0:350) 

DIMENSION  CSTAR(0:350.0:350).CPRIME(0;350)AMS(0:350).BMS(0:350) 

DIMENSION  CMS(0:350)EMS(0;350).FMS(0:350).GMS(0:350) 

2  FORMAT(A) 

WRITE(*.100) 

100  FORMATCINPUT  FILE  NAME  FOR  STORAGE  OF  CONC.  PROFILE:  ',$) 

READ(*  ^)FNAME,QUEST 

OPEN(25JOSTAT=IOSERR=79EILE=FNAME,STATUS=’NEW’) 

79  IF  aOS  .NE.  0)  WR1TE(*  .64)  lOS 
64  FORMAT(‘OPEN  ERROR  ‘,14) 


WRITE(*.%) 

%  FORMATCINPUT  FILE  NAME  FOR  STORAGE  OF  DQDT:  ‘  .S) 
READ(*;Z)FNAME 

OPEN(26JTLE=FNAME3TATUS=’NEWJK)RM=’UNFORMATTED’) 

WRITE(*.98) 

98  FORMATCINPUT  FILE  NAME  FOR  STORAGE  OF  Q/QO:  *.$) 
READ(*^)FNAME 

OPEN(28JTLE=FNAME.STATUS=’NEWJORM=’UNFORMAlTED’) 

WRITE(*.99) 

99  FORMATCINPUT  FILE  NAME  FOR  STORAGE  OF  RO:  ‘S) 
READ(*^)FNAME 

OPEN(30JFILE=FNAME,STATUS=’NEW’J^ORM=’UNFORMATTED’) 

WRITE(*.102) 

102  FORMATCINPUT  FILE  NAME  FOR  STORAGE  OF  QACC:  ‘  .S) 
READ(*^)FNAME 

OPEN(32JTLE=FNAME,STATUS=’NEW’J?ORM=’UNFORMATTED’) 

WRITE(*.104) 

104  FORMATCSHRINKING  DROPLET  RADIUS?  (Y/N):  '.S) 

READ(*.2)  QUEST 

WRITE(*,105) 

105  FORMATCINPUT  KMJCD.KS.THETA  (deg.),  SIGMAAND  LAMBDA:  ‘,S) 
READ(*.*)  KM.KDaCS.THETA,SIGMA.  LAMBDA 

WRITE(*.107) 

107  FORMATCINPUT  JL  (even  #),  10.  Il(even  #).IFIRST  (mulLof  10)’. 

&  /‘ANDILAST(mull.of  10): ‘.S) 

READ(*.*)  JLJO.  I1.IFIRSTJLAST 
WRmB(*.l08) 

108  FORMATCINPUT  DT  AND  IFREQ:  ‘.S) 


READ(*.*)DTJFREQ 
THETAR=2.0*3. 14 1 5926*THETA/360.0 

FTHETA»SIN(THFrAR)*(2.0  +  COS(THETAR))/(1.0  +  C0SCrHETAR))**2 

INC»(ILAST-1FIRST)/10 

DZ.l.(VFLOAr(JL) 

DR»1.0/FLOATaO) 

RO=FLOATai)/FLOATaO) 

R0*1 

ALPHA=(DR/DZ)**2/LAMBDA**2 

WRITE(25.1 10)  KM,KD JCS,THETA.SlGMAa-AMBDA J10.DR.DZ,IFIRST.ILAST.1NC 
no  FORMATCUNSTEADY  STATE  DROPLET  DIFFUSION  IN  A  MEMBRANE.  l.A.D.’, 
&  ‘  METHOD.  WITH  PARAMETERS’// ‘KM  =’.D10.5/‘KD  =  •. 

&  DIOJ  A’KS  =  ‘  J)10.5/.’THETA  =  D10.5.  ‘  deg.  ’/ 

A  ‘SIGMA  = ‘.  D10.5/ ‘LAMBDA  =  ‘  JDIO.S/ ‘RO  = ‘.  D10.5  / 

A  ‘DR  =  ‘.  D10.5/  ‘  DZ  *  ‘  JD10.5/  'IFIRST  =  ‘  J4 
A  /  ‘ILAST  »’.I4/  ‘INC  =  ‘.14) 

DO  10 1=0  J1 
DO10J=0JL 

caj)=o 

CSTARO  J)=  0 
10  CONTINUE 
DO  15  I=0J0 
Ca.0)=l 
CSTARa.0)=l 
15  CONTINUE 
J0=1 

ICOUNT=0 
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N=0 


T=0.0 


DQDT=0 

DQ0DT=0 

DQDTEM=0 

DQDTEIM) 

CM) 

QEEM) 

QEM=0 

QACC=0 

Q0=0 

ROiil 

18  BETA=DR**2/(DT»LAMBDA**2) 
AMS(0)=0.0 

BMS(0)=2.0*BETA  +  4.0 

CMS(0)»-4.0 

DO  20 1=1,11-1 

RI=l 

AMS(I)=1.0/(2.0*RI)-1.0 
BMS(D=2.0*(1.0  +  BETA) 
CMS(I>=-(  1 .0+ 1 .0/(2.0*RO) 

20  CONTINUE 
AMSai)=-2.0 
BMSa  1)=2.0*(  1 .0+BETA) 
CMSai)=0.0 
EMS(0)=0.0 

FMS(0)*6.0*(2.04-BETA) 

GMS(0)=-12.0 


D0  22I=U1-1 


R1=I 


EMS(D=3.0*(1 .0/(2.0*RI)- 1 .0) 

FMS(D=6.0*(1.(>+BETA) 

GMS(I>=-3.0*(l.(kl  .Q/(2.{)*RI)) 

22  CONTINUE 
EMSai)=-6.0 

FMSa  1)=6.0*(  l.Of  BETA) 

GMSai)=0.0 
D023J=1JL-1 
AM(J)=- ALPHA 
BM(J)=2.0*(ALPHA+BETA) 

CM(J)=-ALPHA 

23  CONTINUE 
AM(0)=0.0 

BM(0)=6.0*BETA+8.0*  ALPHA*(1 .0+ICM*DZ) 

CM(0)*2.0*BETA-8.0*  ALPHA 
AM(JL)=2.0*BETA-8.0*  ALPHA 
BM(JL)=6.0*BETA+8.0*ALPHA*(1 .0+KS*DZ) 

CM(JL)=0.0 

IF  (T  .EQ.  0.0)  GO  TO  228 

24  T=:T+DT 

IF  ( T  .GT.  (insert  value))  GO  TO  92 

C  PREVIOUS  LINE  IS  USED  IN  ORDER  TO  STOP  COMPULATION  AFTER 

C  THE  NECESSARY  DATA  HAVE  BEEN  COMPUTED. 

C 

N=N+1 

ICOUNT=ICOUNT  +  I 
Q=Q  +  DQDT*DT/2.0 
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Q0=QOfDQ0DT*DT/2.0 


35 


40 

25 


& 


& 

& 


& 

& 


QED=QED  +  DQDTED*DT/2.0 
QEM=QEM  +  DQDTEM*DT/2.0 
IF  aO  J®.  -1)  THEN 
DO25J=J0JL-l 
DO  35  I=0J1 

DM(I)=ALPHA*CaJ-l)+2.0*(BETA-ALPHA)*Caj)+ALPHA*CaJ+l) 

CONTINUE 

CALL  TRIDAG(0.I1  AMS.BMS.CMS.DM.CPRIME) 

DO  40  I=0J1 
CSTARaj)=CPRIME(I) 

CONTINUE 

CONTINUE 

J=JL 

DO  41  1=0.11 
IF  0  JEQ.  0)  THEN 

DM(0)=-(4 .0<-2.0*  BETA)*CSTAR(0  JL- 1  )44.0»CSTAR(  1 JL- 1  )+2.0* 
(4.0*  ALPHA+BETA)*C(0JL-1M6.0*BETA-8.0*  ALPHA* 
(1.0+KS*DZ))*C(OJL) 

ELSE 

IFO  EQ.  II)  THEN 

DM(1 1)=2.0*CSTAR(1 1  - 1 JL- 1  )-2.0*(  1  .Of  BETA)*CSTAR(1 1 JL- 1 )+ 
(8.0*ALPHA+2.0*BETA)*C(IUL-l)+(6.0*BETA-8.0* 

ALPHA*(1  .Of  KS*DZ))*C(1 1  JL) 

ELSE 

DM(I)=-AMS(I)*CSTAR(I- 1  JL-  1)-BMS(D*CSTAR(I  JL-  l)-CMS(l)* 
CSTARa+ 1  JL- 1)+(8.0*  ALPHAf  2.0*BETA)*Ca  JL- 1  )+ 

(6.0*BETA-8.0*  ALPHA*(  1  .Of  KS*DZ))*C(I  JL) 
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END  IF 


END  IF 

41  CONTINUE 

CALL  TRIDAG(Oai  JEMS34S,GMSDM.CPRIME) 

DO  461=041 
CSTARajL)=CPRIME(D 

46  CONTINUE 
J=0 

DO  45 1=10+141 
IF  a  £Q.  II)  THEN 

DMai)=(6.0*BETA-8.0*ALPHA*(1 .0+KM*DZ))*Cai  .0)+ 

&  (2.0*BETA+8.0*  ALPHA)*Ca  1 .  l)+2.0*CSTARai- 1 .1)- 

&  2.0*(1.0+BETA)*CSTARai,I) 

ELSE 

DM(D=(6.0*BETA-8.0‘ ALPHA-8.0*  ALPHA*KM*DZ)*Ca.O>+ 
&  2.0*BETA*C(I,1)+8.0*ALPHA*C(I.1)-AMS(I)* 

&  CSTAR(M ,  1)-BMS(I)*CSTAR(I.  1  )-CMS(I)*CSTAR(I+ 1.1) 

END  IF 

45  CONTINUE 

DMa0+l)=DM(I0+I)-3.0*(I.0/(2.0*FLOAT(IO+l))-1.0) 

CALL  TRIDAG(IO+14I£MS4^S.GMS,DM.CPRIME) 

DO  47  1=10+141 
CSTAR(I,0)=CPRIMEa) 

47  CONTINUE 
END  IF 

IF  GO  £Q. -1)  THEN 
DO50J=14L-I 
DO  48 1=041 
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DM(I)=ALPHA*Ca  J- 1  )+2.0*(BETA- ALPHArC(U)+ 

&  ALPHA*CaJ+l) 

48  CONTINUE 

CALL  TRIDAG(0J1  AMS,BMS.CMS,DM,CPRIME) 

00491=041 

CSTARaj)=CPRIME(I) 

49  CONTINUE 

50  CONTINUE 
J=0 

00  51 1=041 
IFa^Q.O)THEN 

OM(I)=4.0*CSTAR(1.1)-2.0*(BETA+2.0)*CSTAR(0.1)+ 

&  (8.0*ALPHA+2.0*BETA)*C(0.1)+(6.0*BETA-8.0* 

&  ALPHA*(1.0+KM*OZ)rC(0.0) 

ELSE 

IF  a  .EQ.  II)  THEN 

OMG 1  )=(2.0*BETA+8.0*  ALPH  A)*C(1 1 , 1  )+(6.0*  BETA-8.0* 
&  ALPHA*(1  .Of  KM*OZ))*C(1 1 .0)+2.0*CSTAR(1 1-1,1)- 

&  2.0*(1.0fBETA)*CSTAR(Il.l) 

ELSE 

OM(I)=(6.0*BETA-8.0*ALPHA*(1.0fKM*DZ))*Ca.0)+{8.0* 

&  ALPHA-f2.0*BETA)*C(I,l )- AMS(I)*CSTAR(I-  1.1)- 

&  BMSa)*CSTAR(1. 1  )-CMS(I)*CSTAR(I+ 1 . 1 ) 

ENOIF 

ENOIF 

51  CONTINUE 

CALL  TRIOAG(O.Il4iMS4T^S,GMS4>M.CPRIME) 

00  52 1=0.11 
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CSTAR(U)=CPRIME(I) 

52  CONTINUE 
J=JL 

DO  53  I=0J1 
IFa£Q.O)THEN 

DM(0)=-(4.Of2.0*BETA)*CSTAR(0JL-l)+4.0*CSTAR(lJL-l>+2.0* 
&  (4.0*  ALPHA+BETA)*C(0  JL-  l>f(6.0*BETA-8.0*  ALPHA* 

&  (1.(>+KS*DZ))*C(0JL) 

ELSE 

IF  a  £Q.  11)  THEN 

DM(Il)=2.0*CSTARai-lJL-l)-2.0*(l.(H-BETA)*CSTARaiJL-l)+ 
&  (8.0*  ALPHA+2.0*BETA)*Cai  JL-l)+(6.0*BETA-8.0* 

&  ALPHA•(l.(>^■KS•DZ))*CaUL) 

ELSE 

DM(I>=-AMS(D*CSTAR(I- 1 JL- 1  )-BMS(I)*CSTAR(I  JL- 1  )-CMS(I)* 
&  CSTARa+ 1  JL- 1  >+(8.0*  ALPHA+2.0*BETA)*Ca  JL- 1  )-K6.0* 

&  BETA-8.0*  ALPHA*(l.(>fKS*DZ))*C(IJL) 

END  IF 
ENDEF 

53  CONTINUE 

CALL  TRIDAG(0J1£MSJT^S,GMS,DM.CPRIME) 

DO  54  1=0  J1 

CSTAR(IJL)=CPRIMEa) 

54  CONTINUE 
ENDEF 

DO60I=I0flJl 

RI=I 

DO  55  J=1JL-1 
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IF(  I  .EQ.  II  )  THEN 

DM(J)=2.0*CSTAR(I-U)+2.0*(BETA-1 .0)*CSTAR(U) 

ELSE 

IF  a  £Q.  0)  THEN 

DM(J>:(2.0*BETA^.0)*CSTAR(0J)+4.0*CSTAR(1J) 

ELSE 

DM(J>:(  1.0-1 .0/(2.0*RI))*CSTARa- 1  J)+2.0*(BETA- 1 .0)* 

&  CSTARG  JHl.O  +  1.0  /(2.0»R1))*CSTAR(I+U) 

END  IF 
END  IF 

55  CONTINUE 
J=0 

IF  a. EQ.  II)  THEN 

DM(0)=6.0*CSTAR(1 1  - 1 .0)-«5.0»(BETA- 1 .0)*CSTAR(1 1 .0)+2.0* 

&  CSTARd  1  - 1 . 1  )+2.0*  (BETA- 1 .0)*CSTARa  1.1) 

ELSE 

IF  a  .EQ.  0)  THEN 

DM(0)=I2.0*CSTAR(I,0)4^.0*(BETA-2)*CSTAR(0,0)44.0* 

&  CSTAR(1,1)+2.0*(BETA-2)*CSTAR(0.1) 

ELSE 

DM(0)=3.0*(1 .0- 1 .0/(2.0*RI))*CSTAR(I- 1 ,0)+6.0*(BETA- 1  )• 
&  CSTAR(I.0)+3.0*(l.O+-1.0/(2.0*RI))*CSTAR(I+l,0)+ 

&  (1.0-1.0/(2  •RI))*CSTAR(I-l.l)+2.0*(BETA-1.0)* 

<Sc  CSTARd,  1 )+( 1 .0+- 1 .0/(2.0*  RI))*CSTARd+ 1.1) 

ENDEF 
END  IF 
J=JL 

IF  d  EQ.  0)  THEN 


SO 


DM(JL)=2.0*(3.0*BETA-6.0)*CSTAR(0  JLH  12.0*CSTAR(1  JL)+ 
&  2.0*(BETA-2.0)*CSTAR(0  JL- 1  )+4.0*CSTAR(l  JL- 1 ) 

ELSE 

IF  a. EQ.  II)  THEN 

DM(JL)=6.0*CSTARa  1  - 1  JL)+6.0*  (BETA- 1 .0)*CSTARa  1  JL)+ 

&  2.0*CSTARa  1  - 1 JL-1  )+2.0*(BETA- 1 .0)*CSTAR(1 1  JL- 1 ) 

ELSE 

DM(JL)=3.0*(1 .0- 1  .(y(2.0*RI))*CSTARa- 1  JL)+6.0*(BETA- 1 .0)* 
Sc  CSTAR(I  JL)+3.0*(1  .{y(2.0*RI)+ 1 .0)•CSTAR(I-^l  JL)-(- 

&  (1.0-1  .(y(2.0*RI))*CSTAR(I- 1  JL- 1  )+2.0*  (BETA- 1 .0)* 

&  CSTARG  JL- 1 )+( 1 .0/(2.0*RD+ 1 .0)*CSTARa+ 1  JL- 1 ) 

END  IF 
END  IF 

CALL  TRIDAG(0  JLAM.BM.CMJDM.CPRIME) 

D061  J=0JL 
CaJ>=CPRIME(J) 

61  CONTINUE 
60  CONTINUE 

IFa0£Q.-l)GOTO220 
DO  75  I=0J0 
RI=1 

DO70J=lJL-l 
IF  G  EQ.  0)  THEN 

DM(J)=(2.0*BETA-4.0)*CSTAR(0J)-f4.0*CSTAR(lJ) 

ELSE 

DM(J)=(  1 .0- 1 .0/(2.0*RI))*CSTAR(l- 1  J)-(-2.0*(BETA- 1 .0)* 

&  CSTAR(I  J)-t-(  1.0+1 .0/(2.0*RI))*CSTAR(I+ 1 J) 

END  IF 
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70  CONTINUE 


J=JL 

IFaEQ.0)THEN 

DM(JL)=12.0*CSTAR(1JL)+2.0*(3.0*BETA-6.0)*CSTAR(0JLM.0* 
&  CSTAR(1  JL- 1)+2.0*  (BETA-2.0)*CSTAR(0  JL- 1 ) 

ELSE 

DM(JL)=3.0*(1 .0- 1  .(V(2.0*R0)*CSTARa- 1  JL)+6.0*(BETA- 
&  1.0)*CSTARaJL)+3.0*(1.0/(2.0*RI)+1.0)* 

&  CSTARa+ 1  JL>+(  1 .0- 1 .0/(2.0*Rl))*CSTAR(M  JL- 1 )+ 

&  2.0*(BETA-1.0)*CSTAR(IJL-1M1.0/(2.0*RI)+10)* 

&  CSTAR(I+1JL-1) 

END  IF 

DM(1)=DM(1)  +  ALPHA 

CALL  TRIDAG(1  JL,AM,BM.CM,DM.CPRIME) 

D0  75J=IJL 
CaJ)»CPRIME(J) 

75  CONTINUE 
C 

C  FLOW  THROUGH  THE  MEMBRANE  SURFACE  AT  Z=  1 

C 

220  SUM2=0 
SUM3=0 

D0225I=U1-U 
SUM2=SUM24C(I  JL)*  FLOAT(I) 

225  CONTINUE 

DO  259  I=2JI-2^ 

SUM3=SUM3+C(IJL)*FLOAT(I) 

259  CONTINUE 
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AREA=DR**2/3.0*(C(OJL)*FLOAT(0)44.0*SUM2+2.0*SUM3+C(11JL)* 


&  FLOATOl)) 

DQDT=2.0*3.1415926*KS*LAMBDA*AREA 

Q=Q+DQDT*DT/2.0 

C 

C  ACCUMULATION  WITHIN  MEMBRANE 

C 

SUM2=0 

DO270J=0JL^ 

SUMR2=0 

SUMR3=0 

DO260I=UI-U 

SUMR2=SUMR2+C(IJ)*FLOATa) 

260  CONTINUE 
00  2621=241-2^ 

SUMR3=SUMR3+C(U)*FLOAT(I) 

262  CONTINUE 

AREA=DR**2/3.0*(4.0*SUMR2+2.0*SUMR3+C(IU)*FLOAT(I1)) 

IF(J£Q.O)THEN 

AREAC^AREA 

ELSE 

IF  (  J  £Q.  JL )  THEN 
ARE  AJL=  AREA 
ELSE 

SUM2=SUM2+AREA 
END  IF 
END  IF 


270  CONTINUE 


SUM3=0 


D0  271 
SUMR2^ 

SUMR3*0 

D0264I=lJl-i;2 

SUMR2*SUMR2+CaJ)*FLOAT(I) 

264  CONTINUE 
DO  266 1=2 J 1-2^ 

SUMR3=SUMR3+CaJ)*FLOAT(D 
266  CONTINUE 

AREA=DR**2/3.0*(4.0*SUMR2+2.0*SUMR3+Cai4)*FLOAT(Il)) 
SUM3=SUM3+AREA 
271  CONTINUE 

VOLUME=DZ/3.0*(AREA(>+4.0*SUM3+2.0*SUM2+AREAJL) 
QACC=2.0*3. 1415926*LAMBDA*  VOLUME 
C 

C  EVAPORATION  FROM  MEMBRANE  SURFACE  AT  Z=0 
C 

IF  (KM  .EQ.  0)  GO  TO  320 

M1=I1 

M0=0 

ADD=0 

IF(10£Q. -l)GOTO330 
M0=I0 

ET=FLOAT(I0)/2.0 

IET=ET 

IF  (FLOAT(IET)  .EQ.  ET)  GO  TO  330 
M1=I1-1 
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c 

ADD=DR**2/2.0*(C(M1 .0)+Ca  1 ,0))*(FLOAr(Ml  )+0.5) 

330  SUM2=0 

SUM3=0 

DO  325  I=M04-2>IU 
SUM2=SUM2+Ca- 1 .0)*FLO  ATO- 1 ) 

IF  a  £Q- Ml)  GO  TO  327 
SUM3=SUM3  +  Ca.0)*FLOAT(I) 

325  CONTINUE 

327  AREA=DR**2/3.0*(C(M1 .0)*FLOAT(Ml>+4.0*SUM2+2.0*SUM3+C(M0,0) 

&  *FLOAT(MO)) 

AREA=AREA+ADD 

DQDTEM=2.0*3. 1415926*  KM*LAMBDA*  AREA 
QEM=QEM+DQDTEM*DT/2.0 
C 

320  IFaO  JEQ.-1)GOT0  228 
IF(I0.NE.0)GOTO207 
I(>=-1 
J0=0 
R0=0 

QED=QED+DQDTED*DT/2.0 
DQDTED=0 
GO  TO  228 
C 

C  TOTAL  FLOW  INTO  MEMBRANE 

C 

207  QO=Q+QACC40EM 


C 


C  EVAPORATION  FROM  DROPLET  SURFACE 
C 

DQDTED=KD*FTHETA*3.1415926*LAMBDA*R0**2 

QED«QED+DQDTED*DT/2.0 

C 

C  TOTAL  FLOW  FROM  DROPLET  (THIS  IS  USED  AS  A  MATERIAL  BALANCE  CHECK) 
C 

QT=QOK)ED 

C 

IF  (QUEST  .EQ.  ‘N*)  GO  TO  228 

DUM*3.0/(FTHETA*3. 14 1 5926)*QT*SIGMA 

IF  (DUM  JLE.  1.0)  GO  TO  211 

I0=-1 

J0=0 

R0=0 

DQDTEI>=0 
GO  TO  228 

211  RI(MLO-DUM)**(L0/3.O)/DR 

R0=RI0*DR 
IIOeRIO 

T10=FLOATai0)40.5 

iop=no 

IF  (RIO  .GE.  TIO)  IOP=nOf  1 
IF  OOP  .GE.  10)  GO  TO  228 
I0=I0P 
C 
C 

228  IF  (ICOUNT  .NE.  IFREQ)  GO  TO  24 


5ft 


ICOUNT=0 


WRITE(25, 1  15)DT.T,DQDT,Q.Q0.QEM.QED.QACCJ10 
115  FORMAT(//’DT  = ‘.  D10.5 

&  //’AT  A  TIME  T  =  ‘  J510.5/  34X,’DQ/DT(OLrr)  =  D20.6/ 

&  34X,’Q(OUT)  =‘.D20.6/34X.’Q(1N)  = ‘.020.6  /  34X. 

&  ‘Q(EM)  =‘J)20.6  /  34X.‘Q(ED)  =  *,  D20.6  /  34X, 

&  ‘CKACQ  =‘.D20.6/34X.‘R0  =‘.020.6// 

&  ‘CONCENTRATIONS  ARE’//) 

WRITE(26)T*2220.00.OQOT*44.98 

C  THE  VALUE  THAT  MULTIPLIEO  T  IS  THE  FACTOR  NECESSARY  IN  OROER  TO 
C  CONVERT  OIMENSIONLESS  T  INTO  TIME  IN  MIN.  FT  IS  THE  VALUE  OF 

C  L**2/(O*60).  TO  CONVERT  OQOT  INTO  NG/SQ.CM.MIN  WE  HAO  TO  MULTIPLY 

C  OIMENSIONLESS  OQOT  BY  THE  VALUE  OBTAINEO  FROM  (C0*R0**3*O)/L**2 
C 

WRITE(28)T*2220.00.Q*  100.00 

C  THE  SAME  IS  TRUE  HERE  FOR  T.  BUT  Q  WAS  MULTIPLIEO  BY  THE  VALUE  OBTAINED 

C  FROM  ((C0'*R0**3)/10cm**2)*  10*6  IN  ORDER  TO  OBTAIN  Q  IN  MICRO  GRAM/SQ.CM 
C 

DO80J=0.JL 

WRITE(25.120)  (C(U).  I=IFIRST.  ILAST.  INC) 

120  FORMAT(10(D  10.4 .2X).D  10.4) 

80  CONTINUE 
C 

C  THE  FOLLOWING  WAS  IXDNE  IN  ORDER  TO  CHANGE  THE  TIME  STEP  DURING 
C  COMPUTATION 
C 

IF  ( T  .GT.  0.30 )  THEN 
DT=0.00005 
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IFREQ=800 
GO  TO  18 


ELSE 

IF  ( T  .GT.  0.05000 )  THEN 

DT=0.00001 

IFREQsSOO 

GOTO  18 

END  IF 

END  IF 

GO  TO  24 

92  CLOSE  (25.STATUS=’KEEP*) 

CLOSE  (26.STATUS=*KEEP’) 

CLOSE  (28.STATUS=’KEEP‘) 

CLOSE  (30,STATUS=’KEEP’) 

CLOSE  (32.STATUS=’KEEP’) 

CLOSE  (34,STATUS=’KEEF) 

END 

C 

c 

c 

c 

c 

c 

SUBROUTINE  TRIDAG(IFL-w43.C,D,V) 

DIMENSION  A(0:350)3(0;350).C(0;350),D(0;350),V(0;350) 
DIMENSION  BETA(0;350).GAMMA(0:350) 
BETA(IF)=B(IF) 

GAMMA(IF)=D<IF)/BETA(IF) 
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IFP1=IF+1 


D0  5I=IFP1J. 

BETA(I)=B(I)- A(I)*Ca- 1  )/BETAa- 1 ) 
GAMMA(I)=(D(I)-A(I)*GAMMAa-  l))/BETAa) 

5  CONTINUE 

V(L)=GAMMA(L) 

LAST=L-IF 
DO  10  K=1XAST 
I=L-K 

V(I)=G  AMMA(I)-C(D*  V(I+ 1  )/B  ETA(I) 

10  CONTINUE 

RETURN 
END 

12.2  1*  OlriMnslonal  Model 

c  DIFFUSION  IN  A  PLANE  SHEET  (ONE  DIMENSIONAL  ASPECT) 

C  FINTTE-DIFFERENCE  METHOD  (NON-LINEAR  SOLUTION) 

C  TRIDAG:  SUBROUTINE  FOR  SOLVING  A  TRIDAGONAL  SYSTEM  OF 
C  SIMULTANEOUS  EQUATIONS 

C  CPRIME:  VECTOR  FOR  TEMPORARY  STORAGE  OF  CONCENTRATION 
C  COMPUTED  BY  TRIDAG 

C  DSTAR:  DIFFUSION  COEFFICIENT  AT  OLD  TIME  AS  A  FUNCTION 
C  OF  CONCENTRATION 
C  DDSTAR:  DERIVATIVE  OF  DSTAR 
C  KS  :  DIMENSIONLESS  MASS  TRANSFER  COEFFICIENT 
C  AM.BM.CMJDM  ARE  THE  COEFFICIENT  VECTORS. 

C 

CHARACTER*20  FNAME 
REAL  LAMBDA.KS.M 
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double  precision  ain,bm,cin,dni,c,cprime,(lstar,ddstar 
DIMENSION  AM(0:700).BM(0;700).CM(0;700),DM(0;700). 

&  C(0:700).CPRIME(0:700)J)STAR(0:700)X)DSTAR(0:700) 

2  FORMAT(A) 

WRITE(*.100) 

100  FORMATCINPUT  FILE  NAME  FOR  CONC  FILE:  *.$) 

READ(*^)  FNAME 

OPEN(25JTLE=FNAME.STATUS=’NEW*) 

WRITEr.lOl) 

101  FORMATCINPUT  FILE  NAME  FOR  DATA  FILE:  *,$) 
READ(*;2)FNAME 

OPEN(26fILE=FNAME,STATUS=’OLD’) 

READ(26.*)KS  >1 JL  JFIRSTJL  AST 

READ(26*)DTJFREQ 

WRITE(*.90) 

90  FORMATCINPUT  FILE  NAME  FOR  STORAGE  OF  DQDT:  ‘  ,S) 
READ(*.2)FNAME 

OPEN(27fILE=FNAME,STATUS=’NEW’fORM=’UNFORMATrED’) 

WRITE(*^1) 

91  FORMATCINPUT  FILE  NAME  FOR  STORAGE  OF  Q:  \S) 
READ(*^)FNAME 

OPEN(28JTLE=FNAME.STATUS=’NEW’K)RM=’UNFORMATrED) 

C 

DZ*1.0/FLOAT(JL) 

INC=(JLAST-JFIRST)/10 

WRITE(25.102)KS>I.DTJL45ZJNC 

102  FORMATCKS=  ‘.DIO-SCM*  ‘  J)10.5/’DT=  ‘.DI0.5/ 

&  ‘JL=  M4/’DZ= ‘X>10.5/’INC=  M4//) 
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DO  10J=0JL 
C(J)=0 

DSTAR(J>=EXP(M*C(J)) 

DDSTAR(J)=M*EXP(M*C(J)) 

10  CONTINUE 
C(0)=1 

DSTAR(0)*EXP(M*C(0)) 

DDSTAR(0)=M*EXP(M*C(0)) 

ICOUNT=0 

N=0 

T=0.0 

DQDT=0.0 

0=0.0 

3  LAMBDA=DT/(DZ**2) 

IF  (T£Q.  0.0)  THEN 
WRITE(25.103)  TEAMBDA 

103  FORMAT(‘AT  TIME  T=  ‘.D10.4^X.  'LAMBDA=  ‘,D10.4// 

&  ‘CONCENTRATIONS  ARE’//) 

WRITE(25, 104)(C(J)  J=JFIRSTJLAST.INC) 

104  FORMAT(10(D10.UX)JD10.1) 

END  IF 

4  T=T+DT 

IF  (Q  .GT.  (insert  amount))  GO  TO  92 
C 

C  PREVIOUS  LINE  IS  USED  IN  ORDER  TO  STOP  THE  PROGRAM  AFTER 
C  COMPUTING  DATA  FOR  A  SPECIFIED  Q  OR  TOR... 

C 

N=N+1 


ICOUNT=lCOUNT+l 


AM(1)=0.0 

BM(  1  )=  1 .0+2.0*L  AMBDA*DSTAR(  1 ) 

CM(1)=-(LAMBDA*DSTAR(1)+{LAMBDA/4.0)*(DDSTAR(1))*(C(2)-1.0)) 

D015J=2JL-1 

AM(J)=-LAMBDA*DSTAR(J)+(LAMBDA/4.0)*(DDSTAR(J))*(C(J+1)-C(J-1)) 
BM(J)=  1  ,(H2.0*LAMBDA*  DSTAR(J) 

CM(,>-1  AMBDA*DSTAR(J)-(LAMBDA/4.0)*(DDSTAR(J))*(C(J+ 1  )-C(J- 1 )) 

15  CONTINUE 

AM(JL)=1.0-8.0*LAMBDA*DSTAR(JL)-2.0*DT*DDSTAR(JL)* 

&  ((KS*C(JL))/(DSTAR(JL)*DZ)) 

BM(JL>=3.(>4-8.0*LAMBDA*(DZ*KS+DSTAR(JL))+{2.0*DT*DDSTAR(JL)* 
&  KS*C(JL)*((KS*DZ+DSTAR(JL)-2.0*DZ''KS)/(DSTAR(JL)**2* 

&  DZ))) 

CM(JL)=0.0 
D016J=1JL 
IF(J.EQ.  DTHEN 

DM(1)=C(1}+LAMBDA*DSTAR(1)-(LAMBDA/4.0)*(DDSTAR(1))*(C(2)-1) 

ELSE 

IF(J£Q.JL)  THEN 
DM(JL)=C(n.-  1)+3.0*C(JL) 

ELSE 

DM(J)=C(J) 

END  IF 
END  IF 

16  CONTINUE 

CALL  TOIDAGCl  JL^,BM.CM,DM,CPRIME) 

D017J=1JL 

C(J)=CPRIME(J) 

17  CONTINUE 
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C  FLOW  THROUGH  MEMBRANE  (AT  Z=l). 

C 

DQDT=KS*C(JL) 

Q=Q+DQD'PDT 

C 

DO20J=lJL 

DSTAR(J)=EXP(M*C(J)) 

DDSTAR(J)=M*EXP(M*C(J)) 

20  CONTINUE 
C 

IF  aCOUNT  KE.  IFREQ)  GO  TO  4 
ICOUNT=0 

WRITE(25.105)  DT.TJ-AMBDAJXJDT.Q 

105  FORMAT(//’DT=  ‘X>15.5//’AT  TIME  T=  •X>10.5/34X.’LAMBDA=’ 

&  .D15.5/34X,’DQDT(OUT)=  '.D15.5/34X.’Q(OUT)=  ‘T)15.5 

&  //’CONCENTRATIONS  ARE  7/) 

WRITE(27)T,DQDT 

WRITE(28)T,Q 

WRITE(25.106)(C(J)J=JFIRSTJLAST.INC) 

106  FORMAT(10(D10.4.2X),D10.4) 

C 

C  THE  FOLLOWING  LINES  ARE  USED  IN  ORDER  TO  CHANGE  THE  TIME 
C  DURING  COMPUTATION. 

C 

IF(T  .GT.  10.0)  THEN 
DT=0.0001 


IFREQ=800 


ELSE 


IF  (T  .GT.  0.034)  THEN 

DT=0.0000001 

IFREQ=1 

GOT03 

END  IF 

END  IF 

GO  TO  4 

92  CLOSE(25.STATUS='KEEP’) 
CLOSE(26.STATUS=’KEEP’) 
CLOSE(27,STATUS=’KEEP’) 
CLOSE(28.STATUS=  ’KEEP*) 
END 
C 
C 
C 
C 


SUBROUTINE  TRIDAG(IFXw^.B.C  J).V) 

DOUBLE  PRECSION  A3.C.D,V3ETA.GAMMA 
DIMENSION  A(0;700)3(0:700).C(0;700)JD(0:700).V(0;700) 
DIMENSION  BETA(0:700).GAMMA(0:700) 
BETA(IF)=B(IF) 

GAMMA(IF)=DaF)/BETA(IF) 

IFP1=IF+1 


DOS  I=IFPIJL 

BETA(I)=B(I)- A(I)*C(I- 1  )/BETA(I- 1 ) 
GAMMAa)=(D(I)-A(I)*GAMMA(I- 1  ))/BETA(I) 
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5  CONTINUE 


V(L)=GAMMA(L) 

LAST=L-IF 
DO  11  K=1^ST 
I=L-K 

V(I)=GAMMA(I)-C(l)*Va+l)/BETA(I) 
11  CONTINUE 
RETURN 
END 
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